import math

from pulp import *

points1 = [(2.51,-0.71),
           (1.90,-0.66),
           (2.11,-0.91),
           (2.35,-1.13),
           (2.16,-1.22),
           (-1.34,-1.02),
           (-1.12,-0.55),
           (-0.88,-0.92),
           (1.89,-0.50)]

points2 = [(-3.44,0.87),
           (-3.16,0.72),
           (-3.04,-0.07),
           (-2.40,-1.81),
           (-2.25,-1.69),
           (-1.88,-2.13),
           (-1.14,2.08),
           (-0.94,-2.32),
           (-0.94,2.65),
           (-0.89,2.19),
           (-0.55,-2.26),
           (-0.45,2.59),
           (-0.35,2.39),
           (0.43,1.92),
           (0.46,-2.27),
           (0.70,2.33),
           (0.95,-2.22),
           (0.96,-2.02),
           (1.26,1.43),
           (1.28,-2.00),
           (1.54,1.56),
           (1.55,-1.98),
           (1.57,-1.82),
           (1.61,0.79),
           (1.79,0.05),
           (1.83,-0.96),
           (2.01,-1.29),
           (2.08,1.13),
           (2.11,-0.73),
           (2.41,-0.37)]

points3 = [(-3.44,0.87),
          (3.38,-0.15),
          (-3.34,-0.68),
          (-3.25,1.34),
          (-3.16,0.72),
          (-3.09,2.10),
          (-3.04,-0.07),
          (-3.01,1.40),
          (-2.78,-1.75),
          (-2.40,-1.81),
          (-2.25,-1.69),
          (-2.22,2.34),
          (-2.15,2.17),
          (-2.09,-2.18),
          (-2.03,2.39),
          (-1.88,-2.13),
          (-1.67,2.71),
          (-1.62,-2.19),
          (-1.57,-2.66),
          (-1.52,2.40),
          (-1.37,-2.68),
          (-1.37,2.42),
          (-1.17,2.27),
          (-1.14,2.08),
          (-1.14,2.45),
          (-0.94,-2.32),
          (-0.94,2.65),
          (-0.90,2.46),
          (-0.89,2.19),
          (-0.83,-2.45),
          (-0.74,2.54),
          (-0.55,-2.26),
          (-0.50,2.05),
          (-0.45,2.59),
          (-0.35,2.39),
          (-0.13,-2.60),
          (0.27,2.24),
          (0.43,1.92),
          (0.46,-2.27),
          (0.54,2.00),
          (0.64,2.13),
          (0.70,-2.55),
          (0.70,-2.34),
          (0.70,2.33),
          (0.95,-2.22),
          (0.96,-2.02),
          (1.05,1.72),
          (1.21,-1.86),
          (1.26,1.43),
          (1.28,-2.00),
          (1.32,1.69),
          (1.33,-1.70),
          (1.41,-2.18),
          (1.54,1.56),
          (1.55,-1.98),
          (1.55,-1.58),
          (1.57,-1.82),
          (1.61,0.79),
          (1.79,0.05),
          (1.83,-0.96),
          (1.87,-1.61),
          (1.89,-0.50),
          (1.90,-0.66),
          (2.01,-1.29),
          (2.08,1.13),
          (2.10,0.60),
          (2.11,-0.21),
          (2.11,-0.73),
          (2.11,-0.91),
          (2.16,-0.53),
          (2.16,-1.22),
          (2.19,-0.02),
          (2.28,-0.91),
          (2.35,-1.13),
          (2.39,-0.18),
          (2.41,-0.37),
          (2.51,-0.71)]

points = points1

def distance(p1, p2):
    x1, y1 = p1
    x2, y2 = p2
    return math.sqrt( (x2-x1)**2 + (y2-y1)**2 )

##distances = []
##for p1 in points:
##    distances.append([distance(p1, p2) for p2 in points])

distances = [distance(p1, p2) for p1 in points for p2 in points if p1 != p2]
print 'distances: ', distances

radii = ['r_%s' % i for i in range(len(points))]

indexes = ['%s%s' % (i, j) for i in range(len(points)) for j in range(len(points)) if i != j]
print 'indexes: ', indexes

vars = LpVariable.dicts("x", indexes, 0)
prob = LpProblem("Relaxation", LpMinimize)
prob += lpSum([distances[i] * vars[i] for i in range(len(vars))])


##vars = LpVariable.dicts("Radius", radii, 0)
##
### The LP problem
##prob = LpProblem("Maximum radius", LpMaximize)
##
### The objective function
##prob += lpSum(vars), "Maximum radius sum"
##
### The distance constraints
##i = -1
##for ri in radii:
##    i += 1
##    j  = 0
##    for rj in radii:
##        if j != i:
##            prob += vars[ri] + vars[rj] <= distances[i][j], "Distance %f to %f" % (i, j)
##            j += 1
##
### The problem is solved using PuLP's choice of Solver
##prob.solve()
##
### The status of the solution is printed to the screen
##print "Status:", LpStatus[prob.status]
##
### Each of the variables is printed with it's resolved optimum value
####for v in prob.variables():
####    print v.name, "=", v.varValue
##
### The optimised objective function value is printed to the screen    
##print "Maximum radius sum = ", value(prob.objective)
